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ABSTRACT 

An iteration program is described for the determination of equi- 
librium composition and flame temperature of high energy fuels (HEF) and 
other chemical reactions. 

Tables of thermodynamic functions and equilibrium constants needed 
for the calculations were prepared for sixteen (16) elements and compounds. 
Based on the ideal gas relationship, the method of computing values for the 
tables is shown. The CRC 102A digital computer was programmed to compute 
the values of thermodynamic functions and equilibrium constants presented 


in the tables; thereby providing a reasonably consistent set of data. 


ita 
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]. Introduction. 

The high energy fuels in present use and undergoing experimentation 
produce high temperatures upon combustion. On the assumption that chemi- 
cal equilibrium exists among the products of combustion, it is possible 
to calculate the composition and the flame temperature, thereby pre- 
dicting the theoretical performance of a propulsion system. For a system 
containing N products of reaction the simultaneous solution of at least 
N+] equations is necessary. These may involve (1) mass balance, (2) 
dissociation and (3) energy and entropy balance. For nominal products 
of complete combustion of stoichiometric quantities the solution is com 
plex, but at high temperatures, as dissociation occurs and the number of 
products of combustion increases, the calculations become increasingly 
Gaaricult. 

Since equilibrium composition and reaction temperature are interde- 
pendent, a,Series of successive approximations must be made beginning with 
a trial temperature. A number of methods for solving these equations have 
been proposed, Ref. (1), (8), (9), using trial and error, matrices, or 
other convergent approximations. Hand solutions are laborious, so it was 
proposed to program the solution on the Computer Research Corporation 102-A 
digital computer at the United States Naval Postgraduate School, Monterey, 
California. This project was completed during fall 1958 and spring 1959 
in partial fulfillment of the requirements for the degree of Master of 
Sciences. The problem was based on the constant pressure combustion of a 
boron fuel containing carbon, hydrogen and boron (C,H,B) at atmospheric 
pressure uSing air as oxidizer and including up to 16 products of combustion. 
The program is also adaptable to use with other fuels, explosives or any 


chemical reaction that involves the elements Carbon, Hydrogen, Boron, Oxygen, 





Nitrogen, and Argon. 

Tables containing a consistent set of thermodynamic functions are a 
primary necessity when calculating equilibrium composition and flame 
temperature. For the elements and compounds included in this problem no 
such tables could be found in @ Single source nor in a readily useable 
form. Hence a large portion of the effort in this project was devoted to 
compilation,of data and calculation of the thermodynamic functions given 
in Tables I - XXI. The sensible enthalpy, i.e. enthalpy based on the 
value at absolute zero temperature (H°-H8); the total entropy S°; and the 
change in free energy of formation A Fe , containing enthalpy of formation 
at absolute zero, were tabulated for the temperature range 100° to 9900°R, 
Tables II - XVII. Equilibrium constants were calculated for the 1l reac- 
tions shown and these values tabulated over the temperature range 2000°R 
to 9900°R, Tables XVIII - XXI. All calculations were based on the ideal 
gas relationship. Values of specific heat at constant pressure C° may be 


p 


easily obtained as first differences in enthalpy. 





2. Procedure and Discussion 

As a representative high energy fuel (HEF) containing carbon C, 
hydrogen H, and boron B: triethyl - diborane was chosen for use in the 
calculations. When burned with air containing twenty-one percent oxygen 
(21% 0.), seventy-eight percent nitrogen (78% No) and one percent argon 
(1% A). a number of compounds could be produced. Including products of 
dissociatiqn, 16 elements and compounds, Table I, were chosen to work 
with. Based on methods in Refs. (1), (2) and (3), combinations of these 
products of combuStion produced ll possible equilibrium reactions, with 
some repetition shown by the water-gas reaction. 

A readily available source of tabulated thermodynamic functions did 
not exist so it was decided that a reasonably consistent set of tables of 
data would be produced. Formula mass in Table I represents the latest 
values from the International Standards adopted in 1956. The enthalpy of 
formation at absolute zero temperature (HE Dos Table I, was taken from 
Refs. (4) and (6). A program was written for the CRC -102A digital com- 
puter such that basic values of thermodynamic properties could be stored 
in the computer memory and these used to compute the properties at inter- 
mediate temperatures. The enthalpy function was computed for a series of 


small temperature increments by the following formula from Ref. (1): 
H =H, + C, (T-T, ) 


where the subscript b refers to a value at an arbitrarily selected basic 
temperature near the temperature of concern. Since enthalpy was based on 
the value at absolute zero temperature, the tabulated function is a 
"sensible enthalpy" and is a standard value at any given temperature minus 


the standard value at absolute zero, H°*-H3. 


The entropy was computed in much the same manner by the formula, 


Ref (1) 


ce 


Ss = Sh + C., ln a 


but since the entropy at absolute zero iS equal to zero, the tabulated 
values become the standard value of total entropy S°, at any given ten 
perature. As will be shown later in the flame temperature calculations, 
the fact that entropy equals zero at absolute zero makes this a good ref- 
erence temperature to use during these calculations. 

Basic values for Tables II - XVII were found in various dimensionless 
and dimensional forms in Refs. (4) - (9) which represented the most con- 
sistent sources available. In the computer these values were all con- 
verted to the English Standard of Engineering units (BTU/#mol) and all 
computations made and results presented in theSe units. 

Enthalpy H, entropy S, and specific heat at constant pressure Cy for 
A, CO, CO,, H. Hy. N, Ny O and 0, were taken from Refs. (4) and (5): for 
B. BO, B,03, NO and OH from Ref. (6): and C (B sraphite) from Ref. (7). 
Basic values were taken from the references at temperatures of 100°K, 
298.16°K, 400°K, 500°K, 700°K, 1000°K, 1200°K, 1500°K% 2000°K, 2500°K, 
3000°K, 3500°K, 4000°K, 4500°K and 5000°K. Smaller increments were chosen 


in the low temperature range because that is the region in which C_ changes 


Pp 
more rapidly. A graphical extrapolation was made for Carbon to obtain 


basic values at 4500°K and 5000°K. 


The specific heat C_ may be obtained by taking the first difference 


Pp 
in enthalpies, ieee for CO, Table VII: 


= = = LB 
H = Heeggon ~ Hs4g0en ~ Sp @7 


(44,834 - 43,944) BTU/#mol = C., (5500-5400) 


891 BTU/#mol 
C= 100 °R = 8.9] BTU/(#mol °R) 





The change in free energy of formation was calculated on the basis of 
formation from the elements and includes the enthalpy of formation at 
absolute zero temperature, Table I, hence tabulated values are standard 


values (A Fe ) for any given temperature. From Chapters 7 and 12 Ref. 


(1) 


Fs 


% AF, 


H - TS 


AH- A (TS) 
but since calculations were made on the basis of enthalpy and entropy 
functions at a Single temperature, the term on the right of the equation 


becomes TAS and 

AF = AH® -~ TAS® 
The function AF is the Gibbst function or Gibbs! free energy, named 
after Willard Gibbs who first established this useful relationship. As 
an illustration; for the formation of atomic nitrogen from the free ele- 
mental molecular nitrogen 


1/2 N, =N 


OFe = Z\H° = TAS" 


= (Hy - 1/2 Hy, + (He )- = T(Sy - 1/2 Si) 


Tabulated values of enthalpy, entropy, and change in free energy of form- 
ation are presented in Tables II - XVII. 

The equilibrium reactions in Tables XVIII - XXI were formulated and 
equilibrium constants calculated on the basis of methods found in Refs. 
(1), (2) and (3). An important use of the Gibbs' free energy function be- 


came apparent in these calculations. 


AF. =-RI In K 


which states that when reactants in their standard states are changed to 


4 





products in their standard states there is a change in free energy equal 
to -RT times the natural logarithm of the equilibrium constant, Ref. (2). 


Then o ai 
f° - (AH '- TAS 


Pa aia 
_ Gta) 





Ke 


% 


for convenience in practical problems, the actual values of equilibrium 
constants were tabulated rather than a logarithm function. This limited 
tabulation in the low temperature range, but, with the free energy functions 
and the above formula, the equilibrium constant may be easily determined 

for any temperature. 

Values in Tables II ~- XXI were computed on the digital computer and 
punched on IBM cards which were subsequently used in an IBM-402 Accounting 
Machine to print the tabular format. 

With reasonably consistent values of the thermodynamic functions, the 
problem of equilibrium composition and flame temperature was attacked. 

Since these two parameters are interdependent, a trial temperature must be 
chosen and the composition determined; then a check for energy, entropy 
balance must be made and if the temperature was not correct a new trial 
value must be assumed. This procedure must be continued until the correct 
temperature, hence proper composition, is found. 

For the 16 products of combustion, simultaneous solution of 17 equations 
is necessary. An iterative scheme of solution is used on the digital com- 
puter and this is best shown by an illustrative example. First choose a 
trial temperature, 5400°R, and having values of oe Tables XVIII - XXI, the 
interrelation between combustion products and reactants is established and 


equilibrium components may be determined. Burn the fuel completely with 





stoichiometric quantity of air as oxidizer to determine the amounts of 
nominal (major) products. With triethvyl-diborane the following reaction 


occurs 


(iol om / = 
BH, (C,H + 12 05 + (12) oa Ny + (12) ay A= “620054. INH 


ie 2 


5)3 Z 


% + Bo0, + 44,6 a5 + 0.571 A 


Now Since complete burning rarely occurs in an engineering application 


without an excess of oxygen, include CO, Ho and 05 among the nominal pro- 


ducts. The reaction then takes the form 


78 / " 
BH, (CH), + 12 0, + (12) a No + (12) a A =a CO, + b H,0 


feCoBe Os + Genet veanet aT (COCt og Hy uae ole, 


: 23 2 Ze 


where the letters a-h represent the unknown molar amounts. By a material 


balance equations may be written to equal the number of elements involved. 


a CO, Carbon (C) ; a+f=6 

b H,0 Hydrogen (H) 2b + 2g = 18 

Cc B20, Boron (B) 2c = 2 

d No Nitrogen (N) 2d = 44.6 
e A Argon (A) e = 0.571 
f CO 

e HH, 

h 0, Oxygen (0) Vy ee) 


Since only six (6) equations exist for these eight (8) unknowns one mst 


resort to equilibrium equations for solution of the other two. To solve 





for f choose the following equilibrium reaction from Table XVIII 


SEER, SS 
co, == CO+ 1/2 0, 


then by partial pressure, p, notation 
| y. 
K = Teo + (Po) 
p 


fe 


CO, 


% 
The partial pressure of a gas may be represented by the mole fraction 


times the total pressure. The unknown number of moles of carbon monoxide 
has been assigned the value f; let, sigma, > , equal the total number of 
moles of combustion products, then the mole fraction of CO is f/% . 


With the value of a at 5400°R from Table XVIII the above equation becomes 


ee 
g/s Phe P} 
afs FP 


K. = 0.32898 
p 


0.32898 = 5 (hz P)* 


From the carbon (C) material balance, a = 6 - f and 


2.329 £(tle PY 
6- 


6 (0.329) = 1.974 
h/s P)% 40329 (b/s P)*+0.329 


where P is the total pressure in atmospheres. Solution for g, the moles 
of H.. is completed in a similar manner. Then for each new product of 


combustion a new equilibrium equation must be solved. This points out 





the need for values of equilibrium constants at various temperatures. 
It should be noted that the solution for f contains the unknown 
h, for the number of moles of unburned oxygen; this will also occur in 
the solution for g. Since simultaneous solution of this large number of 
equations is difficult these are best solved by an iterative process. A 
trial value of h is chosen; solution for the other unknowns is made on 
the basis of this trial, then by the summation shown for the oxygen balance 
the trial value is checked. Considerable assistance is gained by making 
a table to show the operations involved, (see Iteration Table page 10). 
The values of _ predict negligible amounts of Boron and solid 
Carbon for this example. For products of dissociation involving other 


minor components, assign the following unknowns. 
j = OH, k =H, 1 =0, m=N, n= NO, q = BO. 


For the first iteration, only major products are considered to exist 

hence all minor products will be zero. The role played by pressure is 
indicated in the Iteration Table, but for convenience was taken as one (1) 
at~osphere. 

At the end of the first iteration the value of h is determined and 
compared with the trial value. The discrepancy in h is taken as the cal- 
culated value minus the trial value. For the second trial value of oxygen 
divide the discrepancy by two (2) and add this to the first trial value 
of h. Now, however, all minor products of combustion should be included. 

After the second iteration extrapolate (or interpolate) the dis- 
crepancy to zero to aid in the estimate of a trial value of h for the third 
iteration. The third trial should then give a reasonable answer provided 


no errors in arithmetic were made. Once programmed, this sort of iterative 





ITERATION TABLE 


h = trial O 
? 


= = 38,87+h+j+k+l+mtn+a 
(n/s p)’* 


” Th/s Py 0.329 


a= 6 ~ f 
a 9 
B” 21.96 (n/s Pee 1 


Perto-- & - J wk 
e=l+q 


fu. 
i 


22.3 = Me- nN 
0.571 (const. ) 
(e/x Py’ 
P 0.0476 b 
(p/5 P)"2 
k =/p (0.0248) (g/s P)” 
1 =2/P (0.0144) (h/s Pye 
2/P (0.00137) (d/s py 


1) 
ii 


nm = 
n = 0.1216 4% H* 
Y 
(c/E. PY* 
= 0.00008 ee ee 
q = 2/P ( 9) (h/s P)# 
2at+b+3c+f+j+l+n+q 


= 
h ( 5 


Discrepancy in h = ealculated 
~ trial 


10 


Za00 


MANES 


» 24,6 
3.440 


2.560 
1.410 


7.590 
IUE 01810. 


22.330 


Up) 


oo f= oo 


Piette 


-0,08 


220 


tele 3S 


244 
32450 


2.550 
1.420 


7.580 
1.000 


Pas} S010 


Ons 71. 
05155 


1.960 


i2oe 
Zi 
0.043 
On 720 


0.001 


2.50 


46,704 


ae Bt 
36530 


Pa FhO) 
1.484 


4. 336 
1.000 


PANN 35) 


0.272 
Oc 17s 


1.160 


1.310 
1.300 
O.083 
0.890 


Q0.001 





process is quickly and accurately completed on a digital computer. 

Once the composition in number of moles of each product of com- 
bustion has been determined, the trial temperature must be checked. 
This is accomplished by an energy, entropy balance. If the temperature 
is incorrect a new value must be chosen and a new composition determined. 
It is advisable to choose the second temperature value such that it will 
be on the opposite side of the true value from the first trial. Then 
a more accurate eStimate may be made and the theoretical flame tempera- 


ture determined. 
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TABLES OF 


THERMODYNAMIC PROPERTIES AND EQUILIBRIUM CONSTANTS 
Symbols 


Gibbs' free energy of formation, BTU/#Mol 
Total heat content, Enthalpy, BTU/#Mol 
Enthalpy of formation, BTU/#Mol 
Equilibrium constant, Dimensionless 
Formula mass, #/#Mol 

Gas constant, 1.9872 BTU/(#Mol °R) 
Entropy, BTU/(#Mol °R) 

Absolute temperature, °R 


Superscript, denotes standard value at given temperature, 
ieee, H° is standard enthalpy 


Subscript, denotes temperature, io@€«, H$ is Standard enthalpy 
at O°R 
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TABLE I 


FORMULA MASS AND ENTHALPY OF FORMATION 


AT ABSOLUTE ZERO TEMPERATURE 


Element or 


Compound M (HE), 
A 39.944 0 
B 10.820 0 
BO 26 «820 - 10,800 
B, 0, 69.640 - 543,600 
@ 12.011 0 
CO 28.011 - 48,963 
CO, 44.011 - 169,143 
H 1.008 92,916 
H, 2.016 0 
H,0 (g) 18.016 - 102,793 
N 14.008 153,216 
Ny 28.016 0 
NO 30.008 38,659 
0 16.000 105,455 
0, 32.000 0 
OH 17.008 18,000 
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